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TECHNICAL MEMORANDUM X- 64581 


ANALYSIS OF A TWO-DIMENSIONAL MASS TRANSPORT PROBLEM 
IN CONTAMINATION OF SPACE FLIGHT EXPERIMENTS 

INTRODUCTION 


The problem of contaminant particles becomes ever more critical as 
the state of the art advances in space flight experiments involving measurement 
of optical and ultraviolet radiation. A method of determining contributions from 
local sources is needed to optimize system performance. Sources such as 
spacecraft dumps and rocket firings are relatively easy to analyze and isolate 
from critical surfaces because they are external and discontinuous in time. 
Outgassing, however, occurs internal to the experiments and continuously in 
time. Quantities of experimental data on outgassing rates are available, but 
there has been no effort to apply boundary conditions to make such outgassing 
data meaningful in a flight situation. The purpose of this study is to make such 
an application and thereby contribute a vital unit to an optimization program for 
design of optical flight experiments. 

It may be worthwhile here to place this problem in context and point 
out the interfaces with other subprograms in the optimization program (Fig. 1) . 
The output of the first subprogram will be source rates. This will serve as 
input to the mass transport subprogram along with sticking coefficients of 


INTERFACE INTERFACE 



Figure 1. Subprograms in the contamination program. 






boundary surfaces. The output of the mass transport analysis will be a mass 
flux at the optical surface of interest. This will serve as input to an optical 
degradation program to determine the effects of this flux upon the experiment. 
All areas shown in Figure 1 are being worked on. 


ANALYTIC APPROACH TO THE TWO-DIMENS I0NAL 
CONTAMINATION PROBLEM 


The eventual goal of this study was to provide a general design tool for 
optimizing optical experiments in space. At the time the study was begun the 
most pressing need for such analysis existed in the Apollo Telescope Mount 
(ATM) program. For this reason the parameters chosen are meant to 
exemplify a typical ATM experiment. Even more critical is the altitude for 
which the calculations are carried out. This altitude determines whether the 
problem is to be purely geometric or is to involve particle -particle collisions 
within the experiment itself. At the ATM altitude of 425 kilometers (235 n. mi) 
the mean free path of a particle is approximately 10 4 meters [1], which 
essentially dictates a geometric approach. Although pressure in the canister 
builds up to higher values the lowest value of mean free path with the doors 
open is approximately 80 meters, which is still a geometric problem. How- 
ever, the doors are closed during the night cycle, since ATM is a solar experi- 
ment, and the pressures build up to values that would result in intermolecular 
collisions. With the doors closed there is no "wind" and consequently no 
momentum transfer; therefore, the problem is still basically geometrical in 
character. The primary effect might be expected to come from the sticking 
coefficient, but this turns out to be surprisingly pressure independent, 1 and 
for this reason the increase in pressure during the night does not seem to 
seriously invalidate the geometric approach. At worst, an effective increase 
in source rate should handle this problem. There is a compensating factor in 
that absorption on the lens is less in the absence of ultraviolet radiation. 


Calculation of the Probability of Particle Striking Optical 
Surface After One Collision with Wall 

To solve the problem analytically it is necessary to obtain the 
probability that a re -emitted particle at point x on the wall will strike the lens 


1. Conversation with J. F. Scannapieco, General Electric, Valley Forge, 
Pennsylvania, July 1969. 
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( Fig. 2) . This is simply the ratio of the angle subtended by the lens to the 
angle into which the particle is emitting. (For three dimensions, it is the 
ratio of solid angles.) This is with the assumption of isotropic emission. 

Such an assumption will be replaced by a more realistic distribution as soon as 
the geometric probability is obtained. 



Figure 2. Basic geometry to be considered. 

General symbols are defined as follows: 

N the number of particles per second from the source at 

( 0 , 0 ). 

N the number of particles per second from the wall at 

x / \ 

(x,y) 

co angle subtended by lens referred to the source (0, 0). 

co' 7 angle subtended by lens referred to point (x, y) 
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X 0 , X 1# Y 0 , Yj corner points of two-dimensional box 
y 0 , y 1 boundary of lens 

1, 2, 3, 4 the number of the wall for later clarity. 

From Figure 2 it is seen that 

x = Y tan 6 , ( i) 

and, from the defined quantities, it follows that the number of particles per 

CO 

second which will strike the lens directly from the source is — N, while 

o, 2ir 

from wall 1 there will be — N particles per second striking the lens. To 

7 r x 

determine some value for N x> it is necessary to scale up the part of the 
drawing concerned with the top wall (Fig. 3) . 



( 0 , 0 ) 


Figure 3. Enlargement of top wall geometry. 


From Figure 3 it can be seen that 


dN = N — 
x 2ir 


9 


which gives the number of particles per second hitting at x in the cone d0 . 
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(Notes 27 t radians is to be replaced by 4 ir steradians in the three-dimensional 
problem. ) The total number of particles that hit between x and x+Ax per 
second is defined by 


N 


x,Ax 


x+Ax 0(x+Ax) 

/ d V^ / ^=2? [•<***>-•<»>] 

X 0(x) 


( 2 ) 


From equation (1) , it follows that 

-l 


0(x) = tan 


(f) 


0(x+Ax) = tan 




A series expansion [2] for the arc tan a is 


-i 111 

tan a = a - — a + —a - — a + . . . ( a < 1) 
o 5 ( 


7T 1 11 1 

2 a + 3 ot 5q' + *”’ QI> ^ 


(3a) 

(3b) 


Using the formulas to first order in Ax for x < y allows the calculation of 

2n 


-■•'(“) '-»■ 


Ax for x < y . ( 4) 


Consider the case x > y 



Here it is necessary to use the binomial theorem, 


n 


lim ( 1+Ax) « 1 + nAx , 

Ax-~0 


to obtain 


tan 




Therefore, to first order in Ax, equation (5) is obtained: 


tf 

h 

z <-« n+1 0 

\2n 

o 



Ln=l V 

V 


Ax for x > y 


(5) 


It is possible to put these expressions in a closed form by observing that 

d 


tan 1 (x+Ax) - tan 1 (x) 

A A Ax 

Ax -*-0 


dx 


-l 

tan (x) 


and 


d -i 1 du , . 

— tan u = - — r — or d tan 
dx i+u* dx 


*(:)]■ Kiri 

\ y 


T~ 
3 


By rewriting equation (4) in a simple form through the use of complex notation 
and comparing this to the geometric series, it can be noted that, for x < y, 


'■"(f) — 1;) 


Ax 

y 


i 

n=0 l 


(r)f 


Ax 1 


1 - 


\2 


yAx 

= -T— S - Atan 
x + y 


If)- 
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and for x > y 



From equation ( 2) , it follows that 


N A = 

x, Ax 


_N 

27T 


# , 

x^ + y 


( 6 ) 


for all values of 



. As mentioned earlier, this number is for assumed 


isotropic re-emission off the wall. This is not a very good assumption and 
will now be replaced by a directionally dependent probability factor. A cosine 
dependence which is illustrated schematically in Figure 4 will be used. 


To make an application of this directional dependence to the top wall 
shown in Figure 2, it is necessary to analyze the geometry in more detail 
(Fig. 5) . The definition of several parameters shown in Figure 5 are 


» = Il+Xo 
y 2 


9 = tan 


-* ( Yli X 

X n -X 


S = lyt - y 0 l cos 9 ; 
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COUNTS 

SEC 



x x + Ax WALL SURFACE 


Figure 4. Illustration of meaning of cos4> dependence. 



Figure 5. Geometry from which the directional dependence is derived. 
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lyi - yo I cos tan 


CO = 


■* ( Ilz y \ 

\Xo-x / 




X„ 


x) 2 + (Yj -y') 2 


f 


(7) 


If it is assumed that all particles are re-emitted from Ax, the number of 

particles per second leaving in the direction 0" is proportional to cos0". 

Denote the number per second leaving Ax in direction 0 by N = N 

9 A x 

cos0 . The sum over all directions must be the total number leaving 

N ;i.e., 
x, Ax 


£ 

2 

/ 

“£ 

2 


N d0 = N A = N 
Ax, 0 x, Ax Ax 


£ 

2 

/ 

-£ 

2 


COS0 d0 = 2N 


Ax 


( 8 ) 


Therefore, 


N 


N 


x, Ax 


Ax 


and N = — N cos0 
Ax,0 2 x, Ax 


is the normalized distribution that gives the angular dependence of particles 
leaving Ax. 

It is necessary to modify the preceding results slightly to take into 
account the fact that some particles permanently stick to the walls and are not 
re-emitted. For this purpose a new parameter P will be defined to be the 
probability that a particle will be re-emitted from the surface. This P may 
be a function of temperature or other variables, but it is not necessary to 
specify such dependence at this point. The particles leaving the surface are 
then reduced by a factor P. It should be noticed that i-P would correspond 
to the percentage (or probability) of particles sticking permanently to the 
surface and will be called the sticking coefficient. Under equilibrium condi- 
tions it is then possible to specify the number of particles emitted from 
Ax(at x) into the angle cu(at 0) . This is the number of particles leaving 
Ax(at x) that will hit the optical surface, 


-r 




p 

2 


N A cos0" 
x, Ax 


(s/r) 

7T 


In two dimensions the number of particles leaving the source, going to x, 
being re -emitted randomly with d' Lambertian distribution and hitting the lens 
is 


N Ax . 

source — lens 


— -N A cos0 
2 x,Ax 


Ivi-Ynl COS 0' 

ir[(X 0 -x) 2 + (Y r y') 2 ] 


Although the cosine dependence is explicit in cos0", it is preferable to 


employ the relationship cos 


(H- 


sin 0 ' and use the expression for 0 ' 


given in equation (7) . The expression for N for the particles hitting in 

the region Ax has already been given [equation (6)] in terms of the source 
rate and geometry. Substituting these expressions into equation (9) yields 



To consider all the particles from the source which make only one collision 
with the "top" (or "bottom") wall and proceed to the lens it is necessary to 
integrate over the top (or bottom with appropriate geometry) wall: 



Before treating the case in which the particle makes more than one collision 
with the walls and then hits the lens it is necessary to consider the "back" 
wall (wall 3) . The geometry for this case is shown in Figure 6. 
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Figure 6. Geometry of back wall contributions. 

The distance from the Ay region on the back wall to the center of the 
lens is given by 

r = \J x’ 2 + (y' - y) 2 = n/ x' 2 + y" 2 , 


and the angle $ is given by 



As previously determined the number of particles per second leaving the Ay 

"area" is PN . , where N . has been defined for the back wall 
y> Ay y, Ay 

analogously to N for the top wall. The number leaving in the direction 

X 3 ZaX 

$ is proportional to cos<3? . After normalizing, as in equation ( 8) , this 
number is seen to be 


11 



cos $ 


N 


Ay,$ 


— ■ N 
2 y,Ay 


For random $ the probability of a particle hitting the lens is 

!y r y 0 l cos tan ^^7 j 


N 


Ay-dens 


= N ) = N 

Ay,$ \ ir ) Ay, 




^<7 x ’ 2 + y " 2 


When the last two equations are combined, the number of particles leaving the 
source, hitting Ay, and re -emitting to strike the lens is seen to be 

p lyi-yol cos 2 tan " 1 

N Ay = — N ■■ J-A- 

source -* lens 2 tt y,Ay ^ x >2 + y M 2 

The number of particles leaving the source and striking Ay is found 
in the same manner that the number striking Ax was determined. Thus, 



N A 

y> Ay 


N x' Ay 

27T x ,;J + y 2 


This is substituted into equation ( 12) . To find the total contribution from the 
back wall of particles making only one collision and proceeding to the lens, 
the sum over all Ay's is taken. In the limit this sum becomes the integral 


K . 1 
back — 


= PNlX„-X 1 lly,-yJ 
lens ( 2 t r ) 2 


Yi cos 2 tan 

; 

Y 0 


V Xn-Xj / y . 

A 2 - Iv'-v) 3 / 2 


[(Xo-Xj)^ - (y'-y) ] 


(13) 


The one -collision contributions from the top and back walls have been 
determined. The bottom wall is simply treated as a top wall and equation (11) 
is used with appropriate geometry. Since only straightline paths are assumed 
there can be no single -collision contribution from the front wall. The con- 
tribution from the source directly to the lens will now be calculated. The 
geometry is shown in Figure 7. 

The direct contribution to the lens is proportional to the angle sub- 
tended by the lens divided by the angle into which the source is isentropically 
emitting, that is, co/2ir , where 
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(a) 


(b) 


Figure 7. Source to lens geometry. 


u> = 




Ay cos 6 
r 



For a source emitting N particles per second into 2i r radians for Ay « x 0 
this expression is reducible to [See equation (6)] : 


N 

source 


lens 



for Ay « x 0 


(14) 


At this point, calculations are possible and have been carried out for typical 
geometries. These will be compared to a more realistic solution which is 
considered in the next subsection. 


Calculation of the Probability of a Particle Striking the Optical 
Element After Two Collisions with Intermediate Surfaces 


The contribution to the optical surface from particles striking a 
surface element and being re -emitted has been found. All particles on this 
surface element were assumed to have come directly from the source; but 
actually some of the particles on any surface element have reached their 
position after at least one collision with another wall. The special case in 
which only one prior collision is assumed will be treated by the following 
procedure: Let the surface element in question be area Ax on the top wall. 
By treating this element as an optical surface the technique developed in the 
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previous subsection can be used to find the contributions to this element from 
particles making one collision with another wall prior to striking the surface. 

It is necessary to make the identification of Ax with I yt~yo I in equations (11) 
and ( 13) and to modify other parameters to fit the format of these equations. 
Symmetry arguments allow all possible wall combinations to be handled with 
only two basic calculational techniques. Both basic configurations are repre- 
sented in Figure 8. 




(bl 


Figure 8. Geometry involved in determining number of particles 
making one collision with the walls before striking x. 

The geometries shown in Figure 8 must be rotated through 90 degrees 
to be put in a format which is obviously suited for application of equations ( 11) 
and ( 13) . The rotated geometries are shown and labeled in Figure 9. 

It should be obvious from inspection of the above figures that the 
words "top" and "bottom" and "front" and "back" are interchangeable. All 
possible combinations of two wall collisions are handled by equations ( 11) 
and ( 13) with conversion of the figures to the appropriate formats. 

For purposes of calculation it is necessary to choose some sample 
geometry. The geometry shown in Figure 10 will be that used for all calcula- 
tions in this report. 
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"LENS" 


Figure 9. Conversion of Figure 8 to required formats. 



Figure 10. Values used in all calculations (all distances in centimeters) . 

Calculations have been performed for the configuration shown in 
Figure 9a with the distances given in Figure 10. It is reasonable to assume 
that the one -collision contributions to Ax depend on the location of Ax and 
this is found to be true. Therefore, Ax was located at seven different points 
ranging completely over the top wall. The contributions from wall 3 to Ax 
on wall 2 were found to vary smoothly with x and assumed values from 
0. 000132 PNAx to 0. 000663 PNAx. A curve was fitted to the points such that 
the contribution found by equation ( 11) can be represented as 


15 



= kj(x) PNAx , 


(15) 


N, 


3-*Ax on 2 


where 

k^x) = 0.000703 -0.0000107X+ 0.000000044x 2 
P = probability that particle will re -emit from wall 3 
N = number of particles per second emitted by source 
Ax = area of surface element on wall 2 

N = particles emitted from wall 3 to Ax on wall 2. 

3 -►Ax on 2 

A very good approximation to equation ( 15) is 


N 


wall 3 -►Ax on 2 


« 0.0004 PNAx 


U6) 


Similar calculations for the contribution to Ax on wall 2 from those particles 
which leave the source, strike wall 4, are re -emitted, and proceed on the Ax 
have been carried out as above. The configuration is that corresponding to 
Figures 8b and 9b and the relevant formula is equation ( 13) with appropriate 
parameters. The results of moving Ax along wall 2 are given on the basis of 
nine points. These nine points yielded contributions to Ax varying smoothly 
from 0.000419 PNAx to 0.000798 PNAx. A curve fit program yielded 

N . . = kj(x) PNAx , 

4— Ax on 2 1 


where the symbols are explained above and the dependence on the location of 
Ax on wall 2 is best described by 

kj(x) = 0.0005605+ 0.0000113x - O.OOOOOOlx 2 

A very good approximation to N ^ ^ 2 is given by 


N 


4-*Ax on 2 


« 0. 0007 PNAx 


(17) 
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In a similar manner, symmetry considerations allow the calculation of 
the contribution to Ax from the front wall. The results found on the basis of 
seven points can be summarized by 


N 


i-^Ax on 2 


« 0.00009 PNAx 


(18) 


By extending consideration to those particles which make 2, 3, 4, or more 
collisions with various walls before landing in the region from x to x+Ax on 
wall 2, it is possible to write an expression for the number of particles 
striking Ax per second: 



x+Ax 


x+Ax 

x+Ax 

N A 

x, Ax 

/ 

X 

dN + P 

X 

I 

X 

dN lx +p2 / 

X 


direct 

one 

collision 

two collisions 


x+Ax 


_n 

. . . + P 

x+Ax 

+ P 8 

/ 

X 

dN + 
3x 

f dN + . . . 

J nx 

X 


three collisions n collisions 


The first integral represents the direct contribution to Ax from the source 
and has already been integrated [equation ( 14)] , and the second integral 
represents the one collision contribution expressed in equations ( 16) , ( 17) , 
and ( 18) . A similar notation is adopted for the remainder of the series after 
integration such that 


N A 

x, Ax 


= N 


27T (x 2 +y z ) 


Ax + Pkj(x)Ax + P 2 k 2 (x) Ax + . . . 


+ P° k (x) Ax+ . . . 
n 


(19) 


The maximum value of P=i for which the series diverges represents the 
physically unrealizable situation in which no particles ever stick to the 
surfaces. All other values of P lead to convergent series, assuming that 
k^(x) are well-behaved. It will also be assumed that the following relation 

holds: 


k ,(x) == k (x) 
n+1 n 


( 20 ) 
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Before an attempt is made to sum this series, the situation in which 
there are holes in the container walls will be considered. Although each hole 
must be separately considered as a possible source of contaminant particles, 
it also serves as a sink for particles. If the area of the sinks is denoted by 
a and the area of the container walls denoted by A, the probability of an 

A 

arbitrary particle re -emitting and striking Ax is decreased by . This 

AH - 3, 

is only approximately exact but the error does not justify the effort necessary 
to find the exact effect which would be a function of the walls on which the 
particles, Ax, and the sinks were located. If this factor is taken into account 
in the definition 



and if this probability is used in equation ( 19) , the equation becomes 


N A = N 
x, Ax Ax 


2v(J + y^) + k n* X * P 


(21) 


The relation given by equation (20) establishes an upper bound to the series. 
If the equality is chosen, the series becomes a geometric series in P' after 
factoring out kj(x) P’ from equation (21) . This maximum contribution to 
Ax at x is given by 


N A 

x. Ax 


NAx 


27r(x*+ y z ) 


ki(x) P' 2 

1 - p. 


( 22 ) 


This value can now be substituted into equation ( 9) to give the contribution to 
the lens from the top (bottom) wall taking all collisions into consideration. 
Thus, equation (11) can be rewritten as follows: 


N 


top— dens 


_ PNiyi-Vnl f 

(2ir) z J 


X n 


Y , 2tt Mx)P 

x 2 +yi 1 - P' 


>2 


-i Y.-v' ~l Yi-v' 

sin tan — 1 ^ cos tan — 1 ^ • dx 

Xn-x X n -x 

n/(X 0 -x ) 2 + (Yi - y f ) 2 


( 23 ) 
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In a similar manner equation (21) can be substituted into equation ( 12) 
for the contribution to the lens from the back wall. All collisions are con- 
sidered. This allows equation ( 13) to be rewritten; 


m - PN iyi-yni 7 1 f Xn - X 1 

back-* lens (2ny~ J [(Xo-X^ + (y'-y) 14 

x 0 

, ark.MP- 1 

1 " P ’ ^ n/(X 0 -x 1 ) 2 + (y’-y) 2 


(24) 


A Simpson's rule numerical integration has been carried out for 
equations (23) and (24) . Some value of k*(x) must be decided upon before 
the integral is evaluated. The values in equations ( 16) , ( 17) , and ( 18) add 
up to a value of 0. 0013 for k t . All walls are considered. Because this value 
is low, the integration will be carried out by using several times this number; 
i.e. , 


kj e* 0. 0065 . 

Although the determination of kj is the weakest link in the analysis, further 
work has shown that the above choice of kj agrees very well with data obtained 
by an entirely different approach. 

The different approach to the problem is to try a Monte Carlo simulation 
to be run on a digital computer and compare the results of the analytic and 
simulation methods. The comparison of results will follow the next section 
on the Monte Carlo approach. 

MONTE CARLO MASS TRANSPORT PROGRAM FOR CONTAMINANT 
PARTICLES ORIGINATING AT A POINT SOURCE 


Because of the extreme difficulty encountered in rigorously determining 

a value of k (x) , the extended source function for particles making n col- 
n 

lisions before striking Ax, it was decided to check the approximation made 

with a different approach to the problem. Although the choice of k (x) turned 

n 
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out to be very good indeed, the Monte Carlo program was found to have a 
number of advantages over the analytic one, as might have been imagined. 

These will be discussed after presentation of the program. 

The approach followed is that of generating particles with the required 
probability distributions and following their motions, keeping track of position 
at every collision with the walls. At each collision a probability that the 
individual particle will stick is generated and this is compared with the sticking 
coefficient. If the particle sticks, a new particle is generated; if not, then the 
particle leaves the wall with the required directional distribution and the 
process is repeated. Generating large numbers of particles builds up statistics 
which provide the desired data. 

Since velocities are of no significance for this part of the problem, only 
directions are generated for each particle. Starting at the source, the space 
inside the container is divided into quadrants (Fig. 11) . 




Figure 11. Quadrants to determine first wall hit. 

For each particle a direction is generated and then tested to determine 
into which quadrant the particle was emitted. Once this determination has 
been made it is simple trigonometry to calculate the position where the 
particle strikes the wall. If the particle is to be re-emitted then a new 
determination of quadrants is necessary if it is to be known which wall the 
particle will strike . Such a calculation is necessary for every collision since 
the angles change as a function of position on the wall. Every time the first 
wall is hit, a check is made to find out if the particle hits the lens. Every time 
a particle hits the lens, this fact is recorded and a new particle is generated 
at the source. 
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The equation and figures for a particle impinging on each wall are 
given in Figures 12 through 15. They start with wall 1 and proceed to the next 
higher number. 



Figure 12. Quadrants for first wall re-emission. 




Figure 13. Quadrants for second wall re-emission. 



Figure 14. Quadrants for third wall re-emission. 
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Figure 15. Quadrants for fourth wall re -emission. 

Every particle which is re-emitted is given a cosine distribution which 
maximizes the probability of coming off perpendicular to the walls. However, 
the chance that it may go either to the right or to the left requires a test to see 
which wall it will hit. This is done by testing on the appropriate a and j3 
given in the above figures. For each wall, then, a new calculation is necessary 
to ascertain where the particle hit. This is done by trigonometry and results 
in the following sets of equations: 

For a particle coming off wall 1 at position y with angle 0 

if 0< oi x = X 0 + (y-Y o )tan0 on wall 2 

if a<0<j3 y = y+(X o -X 1 )/tan0 on wall 3 

if 0>/3 x = X o +(y-Y 1 )tan0 on wall 4 . 

For a particle coming off of wall 2 at position x with angle 0 

if 9<a y = Yj +(Xi-x)tan0 on wall 3 

if a<9<p x = x+( Y o -Yi)/tan0 on wall 4 

if 0>/3 y = Y 1 +(X 0 -x) tan0 on wall 1 . 

For a particle leaving position y on wall 3 with angle 0 

if 0<a x = X 0 +(y-Y o )tan0 on wall 4 

if Q!<0</3 y = y +(X 1 -X o )/tan0 on wall 1 

if 0>jS x = X 1 +(y-Y 1 )tan0 on wall 2 . 
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For a particle leaving position x on wall 4 with angle 0 


0< a 

y = Y 0 + (X o -x)tan0 

on wall i 

cz<6<fi 

x = x +( Yi-Y o )/tan0 

on wall 2 

0>|3 

y = Y 0 + (X t -x)tan0 

on wall 3 . 


The primary task of the computer is simply to keep track of where the 
particle is as it generates new random directions for the particle. The data 
are generated by counting particles emitted versus particles which hit the lens. 
Although this is the only output of the analytic program, it is also possible to 
find the mean number of collisions for a given sticking coefficient and to 
separate the contributions by wall. 

The advantages of the Monte Carlo technique over the analytic method 
are twofold. It is just as easy to consider a source on one of the walls as an 
interior source with the Monte Carlo. Because the source is always at the 
origin of the coordinate system, the program decides from the geometry where 
the source is and handles the situation accordingly. A further advantage is that 
the thermal distribution is more easily handled in the Monte Carlo case. The 
sticking coefficient is temperature dependent and this ability may be necessary 
if there are appreciable thermal changes in the walls as a function of position. 

If in the future it is necessary to incorporate baffles into the design, 
this would be almost impossible with the analytic approach. Although it would 
not be an easy task, it could probably be written into the Monte Carlo program. 


COMPARISON OF RESULTS OF THE MONTE CARLO AND 
ANALYTIC CALCULATIONS AND CONCLUSIONS 


Figures 16 through 19 illustrate the agreement between the two different 
methods of attacking the problem. Because of the presence of the "extended 
source" parameter kj in the analytic program it was necessary to have some 
check on the analytic calculations. Although not as necessary, it was certainly 
as desirable to have an independent check on the Monte Carlo simulation. The 
curves were plotted on a Hewlett-Packard Calculator Plotter, Model 9125A, 
driven by an H-P Calculator, Model 9100A. The analytic curves are exact 
and the Monte Carlo curves are fit with an exponential regression program to 
the individual points obtained from the Monte Carlo computer runs. It is felt 


23 




0.0 0.2 0.5 0.8 1.0 

PROBABILITY OF RE-EMISSION 


Figure 16. Sample calculations of top wall contribution to mass flux at 

optical surface . 
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Figure 17. Sample calculations of back wall contribution to mass flux at 

optical surface. 
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Figure 18. Sample calculations of bottom wall contribution to mass flux 

at optical surface. 
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Figure 19. Sample calculations of those hits going from source directly 

to optical surface. 

that the remarkable agreement between two such dissimilar methods of 
solving the problem is sufficient to allow the conclusion that the Monte Carlo 
method has been proven successful in its application. Such a conclusion is 
almost required before the much more complicated three-dimensional Monte - 
Carlo simulation is attempted. 

It should also be pointed out here that with slight modifications the two- 
dimensional program can be used as an integral part of an experimental pro- 
gram to determine actual sticking coefficients. 

The graphs (Figs. 16 through 19) are plotted such that the ordinate is 
always mass flux incident at the optical surface (times the known source rate) 
and the abscissa is the probability of re-emission from the walls. This 
parameter serves as the independent variable with mass flux dependent upon 
its value. Each wall is considered separately. It was not felt worth the 
effort to sum all the walls and compare the total flux, although this is one of 
the outputs of the Monte Carlo program. It should be mentioned that the 
contribution of those particles going directly from the source to the optical 
surface should be determined by the geometry and should obviously be inde- 
pendent of the probability of re-emission. This is seen to be the case in 
the relevant graph. 
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The computer program for the Monte Carlo calculation is listed in the 
appendix. The language is Fortran IV and the program was run on an IBM 1130. 
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APPENDIX 

COMPUTER PROGRAM FOR THE MONTE CARLO CALCULATIONS 


29 



n n, n 


^EXTENDED PRECISION 

READ (2*170) XONE * YONE i XZERO* YZERO* Y TOP* YBOT 
170 FORMAT ( 6F 10 .4 ) 

C INITIALIZE CONSTANTS 

ALPHi = 0.0 

BETA I * 0,0 

DETERMINE WHICH WALL SOURCE IS ON AND DEMOTE BY Q...... IF THE 

SOURCE IS NOT ON ANY WALL THEN Q = 0 .... AFTER DETERMINATION IS 

MADE THE QUADRANT ANGLES ARE CALCULATED 
IF ( XZERO) 2 , 1 » 2 

1 Q= 1 

ALPHI = QT AN { -XONE » YONE ) 

BETAI = QTAN ( -XONE t YZERO ) 

GO TO 10 

2 IF (YONE >4 *3 *4 

3 Q=? 

ALPHI * QTAN ( YZERO* XONE ) 

BETAI® QTAN(YZERO»XZERO) 

GO TO 10 

4 I F ( XONE ) 6 * 5 *6 

5 Q*3 

ALPHI » QTAN ( XZER0» YZERO) 

BETAI - QTANt X ZERO* YONE) 

GO TO 10 

6 IF (YZERO) 8*7*8 

7 0 ® 4. 

ALPHI = QT AN ( YONE * XZERO ) 

BETAI = QTAN ( YONE *XONE ) 

GO TO 10 

8 Q=0 

AONE® QTAN (YONE* XZERO) 

ATWO® QTAN (YONE* XONE) 

ATHRE = QT AN ( YZERO * XONE ) 

AFOUR® QTAN ( YZERO* XZERO) 

GO TO 10 

10 QSTAR =Q 

XL=XZ ERO-XONE 
YL=YONE-YZERO 

C DO DIFFERENT VALUES OF PSTIK USING K IF DESIRED 
11 DO 96 K®3 *4 

WRITE (3*777) 

777 FORMAT ( 1 PSTIK J L DIRECT TOP BACK BOTTOM*) 

PSTIK=1.2-.1*K 

C INITIALIZE THE COUNTERS 

I SUM=0 
L? 1 

LZERO = 0 
LTWO = 0 
LTHRE = 0 
LFOUR » 0 
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C INITIALIZE THE RANDOM NUMBER GENERATOR 

P « « 1832#K/2 « 0 
T * » 462 l*K/2 * 0 
S = « 677 l*K/2 « 0 
R « .3473#K/2.0 
12 D095 J = 1*10000 

C RUN J-MAX PARTICLES WHERE J-MAX IS END OF J DO-LOOP 

C RE-SET COLLISION COUNT TO ZERO 

6200 1=0 

C RE-SET CO-ORDINATES TO ZERO FOR EACH NEW PARTICLE 

X = 0.0 
Y=0«0 

ALPHA=ALPHI 
RETA=BETA I 
Q=QSTAR 

C CHECK THE COLLISION COUNT TO DETERMINE WHETHER THIS IS NEW 

C PARTICLE WITH RANDOM DISTRIBUTION OR COMING OFF A WALL WITH 

C COSINE DISTRIBUTION (UNLESS SOURCE IS ON WALL.. ...IF SO THE 

C IF TEST IN STATEMENT 16 WILL GIVE COSINE DIST ANYWAY) 

14 IF( 1)99*16*90 
16 I F ( Q ) 99 » 30 *90 

C GENERATE RANDOM DISTRIBUTION ZERO TO TWO PI 

30 CALL RANDN(R) 

ARGI=6.2831853#R 

C DETERMINE WHICH WALL THE PARTICLE HIT 

32 IF (AONE-ARGI ) 50*97 *80 
50 IFtATWO-ARGI 160*97*51 
51 X=YONE /TAN(ARGI) 

GO TO 200 

60 IF ( ATHRE -ARGI ) 70 *97 *61 

61 Y=XONE*TAN( ARGI ) 

GO TO 300 

70 I F ( AFOUR-ARG 1)80*97*71 
71 X= YZERO/TAN ( ARGI ) 

GO TO 400 

80 Y=XZERO*TAN(ARGI ) 

GO TO 100 

C CHECK TO SEE WHETHER PARTICLE HIT THE LENS .... IF IT DID THEN 

C GO TO ORIGIN AND START NEW PARTICLE RUN (INCREMENTING 

C COUNTERS AS THIS IS DONE).*. IF!NOT sGENERATE NEW PROBABILITY 

C THAT PARTICLE WILL GET OFF OF WALL AND COMPARE WITH PSTIK 

100 I F ( Y-YBOT >105*101*101 

101 1 F ( Y-YTOP 1102*102*105 

102 L=L+1 

I SUM = I SUM + I 
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WRITE (3 *666) PST I K * J »L * L2ER0 sLTWO * LTHRE * LFOUR 
666 FORMAT ( 2X ♦ F6* 3 *6 ( 2X » 1 5 ) ) 

I F ( 0) 99 * 103 ♦ 104 

103 LZERO = LZERO+1 
60 TO 95 

104 IF (0-2)99*702*106 
702 LTWO = LTWO+1 

GO TO 95 

106 I F ( Q-3 1 99 *704 * 107 
704 LTHRE * LTHRE+1 

GO TO 95 

107 IF( 0-4)99*706*99 
706 LFOUR = LFOUR+1 

GO TO 95 

105 CALL RANDN(P) 

IF(P-PSTIK)95. 108*108 

C DENOTE WALL HIT AND GENERATE QUADRANT ANGLES 

108 Q= 1 

ALPHA- QT AN ( XL ♦ YONE-Y ) 

BET A* GTAN ( XL » YZERO-Y ) 

GO TO 500 

200 CALL RANDN(P) 

IF(P-PSTIK) 95*201*201 
201 Q=2 

ALPHA * QTANIYL* X-XONE ) 

BETA = QTAN ( YL» X-XZERO ) 

GO TO 500 

300 CALL RANDN(P) 

IF(P-PSTIK) 95*301*301 
301 0*3 

ALPHA* QTAN ( XL.Y-YZERO) 

BETA* QTAN ( XL * Y-YONE ) 

GO TO 500 

400 CALL RANDN(P) 

IF (P-PSTIK) 95*401*401 
401 Q*4 

ALPHA* QTANIYL* XZERO-X ) 

BETA* QTAN ( YL* ( XONE-X ) ) 

GO TO 500 

C UPDATE COLLISION COUNTER 

500 1*1+1 

GO TO 14 
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c 

<50 


85 

86 

91 
C 

C 

92 

93 

94 
110 
111 

399 

120 

112 

113 

210 

221 

220 

222 

223 

310 

331 

320 

332 

333 

410 

441 

420 

442 

443 

97 

65 

99 

64 

95 


GENERATE COSINE DISTRIBUTION ZERO TO PI 
CALL RANDN(S) 

CALL RANDN(T) 

IF(T“S>85*85»86 
ARG* 1*5707963*5 
GO TO 91 

ARG=3. 141592* ( 1.0-T/2.0) 

IF (0-1)99*110,92 

DETERMINE FROM Q WHICH WALL WAS HIT AND CALCULATE THE POSITION 

WHERE PARTICLE HIT 

IF( 0-2)99*210*93 

IF (0-3)99*310*94 

IF( 0-4)99*410*99 

I F ( ALPHA-ARG ) 120*97*111 

X=XZERO+ ( V-YONE ) *TAN ( ARG ) 

GO TO 200 
BETA « ARG 

IF(BETA-ARG) 113*97*112 
Y=Y+XL/TAN( ARG) 

GO TO 300 

X=XZERO+< Y-YZERO) *TAN(ARG> 

GO TO 400 

IF (ALPHA-ARG) 220 *97*221 
Y = YONE+ ( XOME-X ) *T AN ( ARG ) 

GO TO 300 

IF (BETA“ARG)223*97*222 
X=X-YL/TAN( ARG) 

GO TO 400 

Y= YONE+ ( XZERO-X ) *TAN ( ARG ) 

GO TO 100 

IF (ALPHA-ARG >320*97*331 
X=XONE+< Y-YZERO )*T AN (ARG) 


GO TO 400 

IF ( BETA- ARG ) 333 » 9 7 * 332 
Y=Y-XL/TAN( ARG) 

GO TO 100 

X=XONE+ ( Y-YONE ) *T AN ( ARG ) 

GO TO 200 

IF (ALPHA-ARG) 420 *97 *441 
Y=YZERO+ ( XZERO-X ) *T AN ( ARG ) 

GO TO 100 

IF (BET A-ARG 1443 *9 7*442 
X=X+YL/TAN( ARG) 

GO TO 200 

Y=YZERO+(XONE-X)*TAN(ARG) 

GO TO 300 
WR I TE ( 3 » 65 ) 

FORMAT ( ' PARTICLE HAS HIT THE CORNER - 
GO TO 95 
WRITE (3 *64) 

FORMAT ( '.VARIABLE HAS IMPOSSIBLE VALUE 
CONTINUE 


RECYCLE* ' ) 


* ) 


33 



C CALCULATE RATES FOR RUN JUST FINISHED 

RATE = FLOAT(L-l) /FLOAT! J) 

CAVG* FLOAT ( I SUM ) /FLOAT ( L ) 

DIRCT = FLOAT (LZERO)/FLOAT(J) 

WFOUR = FLOAT ( LFOUR) /FLOAT (J ) 

WTWO = FLOAT! LTWO) /FLOAT (J) 

WTHRE = FLOAT (LTHRE) /FLOAT (J) 

C PRINT OUT SUMMARY OF RUN FOR EACH VALUE OF PST IK 

WRITE (3*701 ) J * PST I K»CAVG» RATE »DIRCT»WTWO»WFOUR» WTHRE 
701 FORMAT! ' SUMMARY OF ' * 16 * ' PARTICLE RUN WITH PSTIK EQUAL 1 *F6.3 
1///' AVERAGE NUMBER OF COLLISIONS WITH THE WALLS WAS ' » F7 . 2/ // ' THE 
2 FLUX RATE AT LENS FOR ALL PARTICLES WAS '.F8.5.' TIMES THE SOURCE 
3RATE. ' //' ' THE DIRECT RATE FROM THE SOURCE WASUF8.5//' THE RATE FR 
40M THE TOP WALL WAS '.3X.F8.5/' THE RATE FROM THE BOTTOM WALL WAS* 
5 F8 » 5 / 1 THE RATE FROM THE BACK WALL WAS '.2X.F8.5) 

96 CONTINUE 
98 CONTINUE 
STOP 
END 


SUBROUTINE RANDN(Z) 
DATA 1/1/ 

REAL M 
I F ( I ) 1*2.1 
1 1*0 

M=2. 0**20 
X=566387«0 
A*2 • 0**10+3 «0 

2 x*A*X-IFIX(A*X/M)*M 

Z = X/M 

RETURN 

END 
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